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The physical phase of Causal Dynamical Triangulations (CDT) is known to be described by an 
effective, one-dimensional action in which three-volumes of the underlying foliation of the full CDT 
play a role of the sole degrees of freedom. Here we map this effective description onto a statistical- 
physics model of particles distributed on Id lattice, with site occupation numbers corresponding 
to the three-volumes. We identify the emergence of the quantum de-Sitter universe observed in 
CDT with the condensation transition known from similar statistical models. Our model correctly 
reproduces the shape of the quantum universe and allows us to analytically determine quantum 
corrections to the size of the universe. We also investigate the phase structure of the model and 
show that it reproduces all three phases observed in computer simulations of CDT. In addition, we 
predict that two other phases may exists, depending on the exact form of the discretised effective 
action and boundary conditions. We calculate various quantities such as the distribution of three- 
volumes in our model and discuss how they can be compared with CDT. 
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I. INTRODUCTION 

Causal dynamical triangulations (CDT) [ll-Q is an attempt to construct a non-perturbative theory of quantum 
gravity. Rather than postulating the existence of new degrees of freedom or new physical principles at the Planck 
scale, CDT uses a standard quantum field theory method — path integrals — to sum over space-time geometries 
weighted by the Einstcin-Hilbcrt action. The path integrals are regularised by discretisation of space-time geometry 
into piccc-wisc flat manifolds with temporal foliation. Usually, space-time is divided into discrete spatial slices, each 
having the topology of the three-sphere, which ensures global, proper-time foliation consistent with the Lorentzian 
signature of the metric. Each spatial slice is represented as a triangulation of the three-sphere, made of equilateral 
tetrahedra. The tetrahedra from neighbouring spatial slices are then glued together, thus forming a complicated 
4d manifold, with periodic boundary conditions in time direction. This lattice regularisation provides a suitable 
ultraviolet cut-off and simultaneously reproduces classical general relativity in the infrared limit. 

Although analytic calculations do not seem to be feasible in the full 3-1-1 dimensional CDT, the model can be studied 
by means of computer simulations. After the Wick rotation to the Euclidean signature, the sum over geometries can be 
performed by standard Monte Carlo methods developed earlier for Euclidean quantum gravity 043- recent years, 
it has been shown that this computational approach has a potential to bring many interesting results. In particular, 
the existence of three phases has been observed [1] . These phases have different profiles of the three- volume N3 (t) as 
a function of time (slice index) t. Depending on the values of parameters in the Einstein-Hilbert action, the system is 
either in phase "A", in which N^it) fluctuates randomly from slice to slice, phase "B" in which A^3(i) is localised in a 
single spatial slice, or in phase "C" in which a macroscopic "quantum universe" is formed [9l-[l2j. In this last phase, 
the average value of N3 (t) at each time slice t is well described by the following formula: 
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for Itl > 
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where A^4 = J^t-^^i^) ^^e total (fixed) four- volume of the universe; s is obtained by fitting to the results of 
simulations; the centre of mass is assumed to be at i = 0. The last formula means that the universe produces a 
"droplet" of cos'^(a;) shape, and that this droplet extends as ttsNI^^ in time direction. This shape is equivalent to 
the classical de-Sitter solution. By making a connection with the mini-superspacc model [l3l | it has been concluded 

m Refs. Ii,[iMl,[il 

that, when only the three-volume is concerned, the full CDT model effectively reduces to a Id 
model with three-volumes {7V3(t)} as the sole degrees of freedom, and with the following discrete action: 

^ = ^iz. j;gt) + ^2i^^3 W> (2) 

Here ci , C2 are new coupling constants related to those in the full Einstein-Hilbert action. An important fact is that 
although the action ([2|) completely neglects the internal structure of each spatial slice t, it gives an excellent agreement 
with simulations of the full model. 

In this paper we introduce a statistical-physics model which reproduces the de-Sitter phase of the CDT. Our model 
consist of a certain number of particles which occupy sites of a Id lattice, and microstates (configurations of particles) 
are weighted with the factor . We identify the emergence of the de-Sitter universe with a condensation- like 
transition known from similar statistical models [TBI , [fcj . We show (both analytically and via computer simulations) 
that a symmetrised version of the action ^ reproduces the shape of the macroscopic universe observed in CDT. 
We calculate the width (temporal extension) of this universe and show that quantum corrections make it wider as 
compared to the classical solution. 

Moreover, we show that the effective action describes not only phase C of CDT but also phases A and B, in the 
space of the coupling constants ci, C2. In addition, we suggest that two further phases may exist: "antiferromagnetic" 
phase D in which thin spatial slices of extended three- volume are separated by slices of minimal size, and "correlated 
fiuid" phase E which emerges from phase C for large four- volume A'4 as a result of merging boundaries of the cos^(a;)- 
shaped universe. In all these phases we calculate quantities such as the probability distribution of the three-volume 
or the correlation function for different three- volumes. Lastly, we suggest that by determining analogous quantities 
in CDT it should be possible to test whether the effective action ([2]) is valid in all phases. 

II. MODEL 

In our model, we consider a one-dimensional closed ring of N sites, each of them carrying a positive number of 
particles mi > 1, . . . ,mN > 1. The total number of particles is equal to M. We denote the density of particles by 
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p = M/N. The numbers of sites N and particles M correspond to the numbers of spatial slices and four- volume N4, 
respectively, while the occupation numbers {nii} correspond to three- volumes {N3{t)} of spatial slices in CDT. 

We assume that the probability of a microstate P{mi, . . . , mjv) factorizes into the product of two-point kernels for 
pairs of neighbouring sites, 

P(mi, . . . , uin) = g{mi,m2)g{m2, m3)...g{mN-i,mN)g(mN, mi), (3) 

where 

5(m, n) = exp -ci — C2 . (4) 



The kernel g{m, n) plays the role of a reduced transfer matrix between neighbouring slices of CDT. The above choice 
guarantees that the partition function 
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corresponds to that of CDT with the effective action ^ in the limit of large systems. Our choice is however 
symmetric in n,m as opposed to We shall see later that this symmetry is necessary to reproduce full-CDT 
simulation results. 

Equation ([3]) has the same form as the steady-state probability of a recently introduced non-equilibrium statistical 
physics model of particles hopping between sites of a Id lattice OUBl- A key feature of this model is the condensation 
phenomenon in which a finite fraction of particles becomes localised in a small region of the lattice if the density of 
particles p = M/N exceeds some critical value pc- In particular, in Ref. [l6j | the following two-point function g{m,n) 
has been analysed: 



g{m, n) = K(\m - n\)^J f{m)f{n), (6) 

with two functions K{x), f{m) playing the role of surface stiffness and on-site potential, respectively. This model has 
a rich phase diagram which depends on the choice of K{x) and f(jn). We will briefly discuss some results of Ref. 
[l6j because they are important for the model discussed in this paper. Let us begin with defining the grand-canonical 
partition function 

M {nii} i 



in which the fugacity z is determined from 

z dh\ Zn{z) 




N dz 



(8) 



We note that the left-hand side of Eq. ^ grows monotonously with z. Since phase transitions are related to singular- 
ities of Z{z) — limjv-s-oo Z]\[{z) and its derivatives, we are interested in the behaviour of this function as z approaches 
the radius of convergence Zc of Z{z). If Zc is infinite, there is always some z > which obeys Eq. for any p. This 
means that Z{z) does not have a singularity for < 2; < 00 which is the physically relevant range of z. Also, both 
ensembles, the canonical and the grand-canonical one, are equivalent in the thermodynamic limit in this case. The 
partition function Zn{z) can be expressed as 

Z N {^^^ — ^ ^ T^ra-im2^ra^m^ ' ' ''^mjv^^ii Tr2^(z) , (9) 
nil 

where T(z) is a square M x M matrix defined as 

T;n„(z) = z("+")/2g(m,n). (10) 
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If we now define 4>m{z) to be the normalised eigenvector of Tmn{z) to the largest eigenvalue Aniax(-z), 

'^T^ri{z)<f)n{z) = Aniax(-z)'/'m(-g), (H) 
n 

we obtain for large N that Z7v(z) = Aniax(^)^- We ean also calculate the probability p{rn) that a randomly chosen 
site has m particles: 

= lim V T,nm.2Tm.^m3---T^,^ni^4L{z), (12) 

The eigenvector (j)m{z) decays with m, and so does p{m). This is guaranteed by the fact that p from Eq. (|5]) is finite. 
In this case, the system has a finite number of particles at every site ~ we say that the system is in the "fluid" phase. 
One can also show that there are only local correlations between different m,;'s in this phase. We shall therefore call 
this phase a weakly-correlated fluid. 

On the other hand, if Z{z) = Iivun^oo Zm{z) has some finite radius of convergence Zc < oo, the derivative in ([5]) 
can either grow to infinity for z ^ Zc, or tend to a finite constant. In the first case, we again have no phase transition, 
because for any p there exists some real z < Zc which obeys Eq. ([5]). However, if Zc < oo and dZ{z) / dz\z^zc ^ const, 
there is a critical density 

Pc = y^TO0r«(gc)^, (13) 

m 

above which the grand-canonical ensemble does not exist. This indicates a phase transition from the fluid to the 
condensed state. 

The nature of the condensate depends on K{x) and f{m) which define g{m^ n). If K{x) — const and f{m) falls off 
sufficiently fast, the condensate spontaneously forms on one randomly chosen site, breaking translational invariance. 



This is precisely the balls-in-boxes (B-in-B) [17[ or zero-range process (ZRP) [18| condensation. We shall note here 
that the B-in-B model has already been successfully applied to the transition between crumpled and elongated phase 
in Euclidean quantum gravity models [l9l . [20l | . 

If K{x) decays with x, the condensate extends to more than one site. The width W of the condensate grows as 
some power a of its volume, W M". The condensate can be either bell-shaped, or rectangular, depending on exact 
forms of K{x) and f{m). We see that this closely resembles the features of the macroscopic universe from phase C in 
CDT. This type of phase, which we shall call a "droplet" phase, will be discussed extensively in Section HVl However, 
if the extension W of the condensate becomes comparable to the linear extension of the system N, both ends of the 
condensate merge and the particles spread uniformly in the system. This phase differs from the weakly-correlated 
fluid phase which exists for p < pc in that the occupation numbers {mi} are correlated. We shall call this phase, 
rather obviously, a "correlated fluid" phase. 

It is important to note that the existence of these phases does not depend on the details of K{x) and f{m) (O, 
which often only affect the shape of the condensate and the value of the critical density. In what follows we shall 
use the analogy between this model and the effective Id model of CDT to study the emergence of the bell-shaped 
quantum universe. We shall also investigate the phase diagram of the model, assuming that the effective action is 
valid in all phases. A small difference between our model and the model from Refs. [ill, 13 is that the two-point 
kernel g{m,n) of the Id effective CDT model (as given in Eq. (U)) has a slightly different form that Eq. because 
K{x) depends not only on the difference between two consecutive occupation numbers m, n but also on their absolute 
magnitudes m, n: 



g{m, n) = K{\m - 7i|/Vm + n)^f{m)f{n) (14) 

However, as we shall see, the only new result is the existence of the antiferromagnetic phase which is not observed in 
the model with kernels of the form (151). 



III. PHASE DIAGRAM 



We begin with presenting the results of Monte Carlo simulations of our model (see Appendix for details), which 
reveal its rich phase structure. Anticipating the existence of condensed/fluid phase, and also extended/localised 
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FIG. 1. Phase diagram determined from Monte Carlo simulations for A'^ = 80 and M — 18100. 



condensates, we define the following quantities which allow us to detect phase transitions: 



E 



M2 



7 = 1- 
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2M 



(15) 
(16) 
(17) 



These quantities assume values between and 1 and play the role of order parameters in the limit N, Af — > oo. The 
parameter a is the inverse participation ratio for site occupation numbers and it measures the degree of localisation: 
for a delocalised microstate in which all to^'s are roughly the same, a « 1/N — >■ for — > oo. However, if one rrii 
is much larger than others, cr — > 1. The value of the parameter 7 indicates whether there are any sites with minimal 
number of particles to^ = 1 in typical configurations (slices with the smallest possible three- volume in CDT): 7 = if 
there are such sites, whereas 7 > if all sites are occupied by larger numbers of particles. The parameter S is related 
to the surface roughness or stiffness of typical configurations and is close to zero for configurations in which {rrii} 
dramatically change from site to site, and close to one for relatively smooth configurations. 

We have used the parameters a, 7 and 6 to determine the phase diagram shown in Fig. [1] (see also Fig.[2]for examples 
of plots of (T, 7, S) in the phase plane of the parameters ci and C2, by simulating our model for fixed N = 80, M = 18100, 
and different pairs of ci, C2. Snapshots of typical configurations in each phase are shown in Fig. [3] Our phase diagram 
includes both positive and negative ci,C2. One might be worried that negative coupling constants should not have 
any physical meaning in the CDT, because the effective action Seg would be unbounded from below for negative ci, C2 
and hence the partition function was ill-defined. However, as we consider here the system with a finite number of 
sites TV and particles M, the action is bounded and the partition function is well defined. 

Looking at Fig. [1] we can distinguish five different phases in the (ci, C2) plane for fixed N, M: 

• Droplet phase: a finite fraction of particles (typically almost all particles) form a bell-shaped condensate extended 
over W ^ 1 sites of the lattice. The shape of the condensate can be approximated by Eq. The droplet 
phase is observed for ci > and C2 > C2,crit(ci), where the shape of the critical curve C2,crit depends also on 
A'' and M. This phase corresponds to the macroscopic universe phase "C" in CDT. The width W and other 
properties of the condensate will be discussed in Section |IVl The values of the order parameters are as follows: 
a is of order 1/W, 7 = and (5 > 0. 

• Correlated fluid: particles are distributed approximately uniformly over all sites of the lattice. The occupation 
numbers fluctuate around the average value {rrii) = P, but the typical size of fluctuations is small as compared 
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FIG. 2. Order parameters: (a) (t(c2) for ci — —0.8, (b) 7(02) for ci = 0.5 and (c) S(ci) for C2 = —0.5. 
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FIG. 3. Typical configurations for all phases from Fig. [T] (a)-droplet, (b)-correlated fluid, (c)-antiferromagnetic, (d)-localised, 
(e)-uncorrelated fluid. 



to the average. This phase is observed for ci > and C2 < C2,crit(ci). In the thermodynamic hmit, we expect 
the order parameters to be cr = (of order for finite system), 7 > and S > 0. 

• Antifcrromagnctic fluid: typical configurations contain alternated occupied/empty (i.e., containing only one 
particle) sites. This phase is observed when both ci and C2 are negative. The number of empty sites increases 
when ci or C2 grow. In the thermodynamic limit, the order parameters in this phase are tr = (of order 2/N 
for finite A^), 7 = and 5 = 0. 

• Localised phase: in a typical configuration, almost all particles occupy a single site, while the remaining sites 
have only small numbers of particles of order 0(1). This phase is observed for ci < and C2 > 0. The order 
parameters are ct = 1, 7 = and (5 > 0. This phase may correspond to phase "B" in CDT. 

• Uncorrelated fluid: Particle occupation numbers are uncorrelated and there is no condensation regardless of 
the density of particles p. This phase is observed in a small region close to the origin of the (ci,C2) plane: 
Ci « 0, C2 « and it may correspond to "A" of the CDT model. 

Interestingly, as we have already mentioned, there are two new phases: the correlated-fluid phase and the antiferromagnetic- 
fluid phase, which have not been observed in computer simulations of CDT. In next sections we shall present some 
arguments supporting the existence of these new phases in the full CDT quantum gravity model. 

We shall now give a crude mean-field argument supporting our phase diagram, based on estimating the value of 
the action 



^ V (mi + m,+i)/2 J 

for typical configurations in different phases, and assuming that, for given ci and C2, the phase with the least value 
of the action is selected. Although we neglect quantum fluctuations of m^'s in this section, we shall see that our 
approach reproduces the phase diagram quite well. Quantum fluctuations will be discussed in the next section. 
The mean-field action for the droplet of width W shown in Fig. [3^ can be approximated as 

where we have assumed that the average shape of the droplet is m.i = {M/W)h{i/W) and that fluctuations can be 
neglected in the limit of large M. We assume that h{x) is fixed and the only degree of freedom is the width W of the 
droplet. Equation HH) can be further simplified if h{{i + 1)/W) = h{i/W) + h'{i/W)/W, 

S^roplet^C,^ dx^^+C2W'/'M'/' l\xh{xy/' . (20) 



(mj+i - rrij)^ 1/3 
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FIG. 4. Phase diagram obtained by comparing the action of typical configurations in different phases. 



The integrals over dx will be explicitly calculated later, now we just treat them as two unknown constants. Searching 
for W which minimizes the action we obtain W - {ci/ C2fl^M^/'^ and, finally, 

5droplet OC c\'^cl'^M^/^ . (21) 

We see that the above calculation predicts the extension W of the condensate to grow as ~ 

Ml/4, gj^^^ 

come 

back to that later. Now, let us consider the energy of the correlated fluid phase (see Fig. [SJj): 

5corr.fluid ~ Nc, ^^-^ ^ + ca^Vpi/s . 22 

p 

Assuming that ~ p + Ami where Ami is of order due to stochastic fluctuations, we obtain 

^corr.fluid (X N Cl + C^N^/Hd^'^ . (23) 

The value of the action for a typical configuration in the antifcrromagnctic phase (see Fig. |3]:) is 

^antifcrr. « 2ciM + C2n^/^M^'^ (24) 

where we have assumed that there are n peaks of height ~ M/n, separated by empty sites. We can use the last 
formula also to estimate the action in the localised phase (see Fig. [3jl) by setting n = 1: 

^localised «2ciAf + C2Af . (25) 

Comparing the values of the action for different ci, C2 and taking the least one, we obtain for large N, M the phase 
diagram shown in Fig. |4l The diagram agrees qualitatively with the experimentally obtained one in Fig. [TJ The 
lines separating different phases are at ci = and C2 = 0, except for a line between the droplet and the correlated 
fluid phase, which has a more complicated shape and will be discussed in Sec. [V] The reader may wonder why we did 
not estimate the action in the uncorrelated fluid phase. The reason is that this phase is dominated by fluctuations 
(entropy) rather than by the action (energy) (|18p which vanishes for ci = C2 = 0. Although this phase exists only at 
a single point (ci,C2) = (0,0) in the phase space in the thermodynamic limit, we expect that for finite systems we 
discuss here, the uncorrelated-fluid phase extends to a small region around (ci,C2) = (0,0). 

We conclude this section with a technical remark. Because our model is motivated by the CDT model of quantum 
gravity, we prefer to use the language of quantum physics rather than that of statistical physics in the paper. If 
one used statistical physics language instead, one would replace the action S by (iE, where /? = ci would be the 
inverse temperature, E would be the energy of configurations, and C2/C1 would be the second parameter (besides j3) 
of our model. The partition function could then be written as Z = e^^^ ~ S{m} i where F would correspond 
to the free energy of the system, including the entropic contribution coming from the sum over all microstates. In 
quantum physics, F is rather referred to as the effective action and the entropic contribution as to the contribution 
from quantum fluctuations. In the next section we shall estimate the contribution from quantum fluctuations to the 
droplet phase and show that these fluctuations lead to the widening of the effective universe as compared to the 
classical de-Sitter solution. 
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FIG. 5. Shape of the universe in the droplet phase averaged over n = 1000000 configurations for Ci = 1 and C2 = 5 and A'^ = 256 
and M — 50000 (blue points), compared to the cos^ shape (Eqs. (|30p and (|3ip ) with the width parameter given by the classical 
solution Eq. p2p W = Wo{M) = 55.5 (red curve), and a more accurate result including quantum and finite size corrections 
Eq. dSIl) = W2{M) = 59.64 (black curve). 



IV. DROPLET PHASE - THE MACROSCOPIC UNIVERSE 



In the droplet phase, which exists for positive coupling constants ci, C2, the condensate takes the form of an extended 
"droplet". In Fig. [5] we show the average shape of this droplet obtained in numerical simulations (see the appendix 
for details). The envelope of the droplet has a cos'^ form and its extension scales as ^ M'^/^ (see Fig. |6]) as determined 
already in the previous section. We will now find the function h{x) and calculate the integrals from Eq. (j20p to find 
the coefficient in the power law W ^ M^/*. Let us first assume that in the limit of large system sizes TV, M — )■ oo and 
p = const, fluctuations of {mi} can be neglected, so that 



rrii 



(26) 



where rrii denotes the average occupation number at site i. The shape of the condensate can be obtained by minimising 
the action 



N 



5({m,}) = ^ 



ci 



(m,+i - TOj)^ 
{fhi + mi+i)/2 



C2m 



1/3 



(27) 



Going into the continuous limit: nii — > m(t) and mi^i — rrii ~^ fTi^t), with m{t) defined on the interval < t < N , 
we see that the following functional has to be minimised with respect to m[t): 



S[m{t)] 



N 



m(t) 



dt. 



(28) 



with an additional constraint that jj^ m{t)dt ~ M . Using the method of Lagrange multipliers we obtain the following 
Euler-Lagrange differential equation for m{t): 



-jrn{t) 



m'{t) 
m{t) 



2ci ^ - a = 0. 

m(t 



(29) 



where a is the Lagrange multiplier used to fix the total number of particles M . This equation is exactly soluble: 

M 



m{t) 



W 



h{t/W), 



(30) 



where h{x) is the "cos^" shape of the droplet, 



h{x) 



^cos3(7r(a;- 1/2)) 



0, 



^sm^{TTx), 0<x <1 

X < OT X > 1 



(31) 
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FIG. 6. The width W versus the number of particles M for different pairs (ci,C2) = (1, 2), (2, 2), (1, 5), (2, 5) (from left to 
right). Black symbols correspond to numerical data, red line shows the classical solution W = Wo{M) (|32|l . black dashed line 
corresponds to the solution W = W\{M) of Eq. (|50|l taking into account quantum corrections and black solid line shows the 
quantum solution including interface effects W = W2{M) from Eq. (|5ip . 



and W is the width of the droplet, 



= Wo(M) = (^^^ « 6.66432 X Ml/4 (^^^ ^32) 

Equations ([5(1)) and ([51]) are equivalent to Eq. ([T} up the position of the centre of mass which is shifted from x = 
to X = 1/2 (we have used the freedom of shifting the droplet to x = 1/2 for the future convenience). The width W 
is uniquely determined by M, ci,C2 and it grows as expected as ~ M^/^ for large systems. Equation ([50]) shows that 
the average height of the droplet M/W scales as ^ M'^/^. Remembering that M plays the role of four- volume of the 
corresponding CDT model, we see that the height is proportional to the three-volume of spatial slices. This is one of 
the reasons why the droplet is considered to be a manifestation of a macroscopic universe in Refs. [9l-[T2|. 

The shape observed numerically closely follows the classical solution (PT|) . see the red curve in Fig. [51 However, the 
width of the droplet W observed in numerical simulations is larger than the one calculated from Eq. ([32[) . as shown 
in Fig. [51 (red curves). The reason is that calculation that lead to Eq. ([5^ neglect quantum fluctuations. We will 
now calculate quantum corrections to W assuming that they leave the shape of the droplet intact. This assumption, 
as we have mentioned, is corroborated by simulations. Our reasoning follows in part the lines of Ref. [l6j in which 
the spatial extension of the condensate has been calculated analytically by splitting the system into two parts: the 
condensate and the fluid background. Proceeding in a similar way, we assume that the total free energy F{W) of the 
system having a condensate of width W can be approximated by 

F{W) « i^backgrou„d(Af ~W) + Fd,.opict(Vt^) = In Z,,,t{N - W, p,{N - W)) + In Zdropiot(T^, M) 

= -W\n\ 

max ■oplct 

{W,M) + 0{N), (33) 

where Zcrit is the canonical partition function for the system with N — W sites being at the critical density, and .^droplet 
is the partition function for the condensate extended over W sites and containing M = (1 — [l — W/N]pcnt/p)M 
particles. If the density is high (the case relevant for CDT), Pcrit/p ^ 1 and we can assume M k, M . Equation 
([33ll states that the free energy of the system is the sum of free energies of the fluid and the droplet, and neglects 
contributions from the boundaries between these two coexisting states. The partition function for the bulk reads 
Zcnt{N — W) ~ A^jjx^, where Amax is the maximal eigenvalue of the matrix Tmn{zc) defined in Eq. ([TU|. The 
partition function of the condensate .^droplet (W^) reads: 



CJO CJO 



^dropiet(M^, M)=Y,---Y1 [Swiir^M})] <5 h\/ - ^ m, , (34) 

mi — Imw— 1 \ i / 

where mo = mw+i — 1 and 

is the action for the droplet of size W . The standard way of estimating the contribution of fluctuations is to expand 
each rrii around its average value rhi, rrii = rhi + Ami , and to assume that the fluctuations Am; arc Gaussian. In this 
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approximation, 



exp 



where Arrii are now continuous variables and the matrix A is the matrix of second derivatives (the Hessian), 

d'^Sw 



A 



'■^ drrtidrnj 



(36) 



(37) 



calculated for {nii} which correspond to the classical solution = M/Wh{i/W) (see Eqs. ([50]) and ^^)- Using the 
integral representation of the Dirac delta 



d{k) 



dq 
2^' 



we obtain 



^droplet (VF) = e 



-Sw[{m}] ^ 
, 27r 



cxp 



— ArriiAijAmj + iq \ 



(38) 



(39) 



We can now calculate the Gaussian integral over Am^'s using the standard result: 



d n cxp 



(2^ 
det A 



cxp 



Wb,b,{A-% 



(40) 



where A ^ denotes the inverse of A. Taking bj = iq for all j we have 



^dropict(W^) = e-^"-[^'""->l 



(2^ 
dct A 



W /■ + 00 



d(7 
2^ 



exp 



and, performing the last Gaussian integral over q, we obtain that 

Fdropiet(M^) - lnZd™piet(M^) = -Sw[{fhi}] + Q(H^), 
where Q{W) correspond to a quantum correction to the free energy: 

Q{W) = ln(27r) - \ Indct A - \ In j 



(41) 



(42) 



(43) 



The first term in Eq. (|42)) is just the action ((35|) calculated along the classical trajectory and it can be easily evaluated 
in the continuous approximation: 



2M/2 



r2/3 



(44) 



where we have inserted m{t) from Eqs. ([50)1 and (|3ip . The quantum contribution Q(M^) to the effective action 
from Eq. ([321) consists of three terms. The first term -^^^ ln(27r) is trivial. The second term — ^Indet^ is more 
complicated because it contains the determinant of A. To evaluate this determinant, we first observe that the matrix 
A is tridiagonal, with only non-zero elements being 

2c2 4ci(-mi_i + TOi)2 8ci(-mj_i + mi) 



Aj. 



4ci 



9^5/3 (mj_i-|-TOi)3 {mi-i+miY rrii-i + 

4ci(-mi + mi+i)2 8ci(-mi + 7T!,,;+i) 



4ci 



4ciVr 



(m^ + mi+i)3 
4ci(-mi + mi±i)2 



(mi + mi+i)2 rui + rrii+i Mh{i/W) ' 
4ci 2cil^ 



(mi-|-mi±i)3 mi + m 



Af/i(Vl^)' 



(45) 
(46) 



11 



We see that Ai,i±i so the determinant det^ can be approximated by Aei A « (deta)]^J^j^ An, where the 

matrix a is a tridiagonal matrix with diagonal elements an = 1 and off-diagonal ones ai^i±i = —1/2. One should 
note that due to the periodic boundary conditions also the corner elements aiN and clni of this matrix should be in 
principle equal to —1/2. In this case the matrix a would have a zero mode. The zero mode has been however removed 
by fixing the position of the centre of mass to be at N/2. With this choice one can safely set oiat = ajvi = leaving 
only the tridiagonal structure of the matrix a. The determinant of this matrix det a = [N + 1)2^^ is independent of 
W , hence the whole dependence of quantum corrections on W is in the factor Hi^i ^a- We can now estimate that 

- 4,c^W Ac^W 128ciW^ 

This is the leading term in Q{W). We shall now argue that the last term J2i ^he quantum correction 

Q{W) can be neglected. The reason is that because Aij oc W/{Mh{i/W)), elements of the inverse matrix A~^ have 
to be proportional to a product of different powers of W, M . Therefore, the sum '^i^^ ^'^so be proportional 

to a certain power of M times a certain power of W (one can show using the fact that Aij is a Laplacian matrix 
times a diagonal matrix with elements ^ l/h{i/W) that Y^j{^^^)ij — MW'^OiVj), and its logarithm will give only 
a sub-leading correction ^ InPF to Q(W^), whose leading behaviour is ~ W^lnW^. 
In summary, the quantum correction approximately reads 

g(I^)^— In- Iniy + lnM +0(A^), (48) 

2 V 64ci / 

and, inserting Eqs. (05]) and (jH]) into Eq. (|^^. and then Eq. (02) into Eq. p3p we obtain the final expression for the 
free energy of the system: 

The width W of the droplet is determined by the maximum of F{W\ Taking a derivative with respect to W we 
finally arrive at an equation for the spatial extension W : 

- In A_ + c,9.^— - c,^-,^ + - fin — + In M - In 1^ - 1 j = 0. (50) 



In the limit of large M, this equation leads to the same expression as Eq. p2[) . For finite M we solve it numerically 
for W . The solution gives a function W = Wi{M) which includes quantum corrections. The maximal eigenvalue Amax 
of the matrix Tm„ from Eq. pop which is necessary to solve Eq. (|50p can be determined by numerical diagonalisation 
of Tmn truncated at to, n w 50. In Fig. Elwe compare W = Wi{M) calculated as a root of Eq. ([50)) and W = Wo{M) 
obtained from the classical formula (j32p . In the same plot we also show values of W measured in simulations of 
the model for different ci,C2. We see that the solution W = Wi{M) which takes into account quantum corrections 
reproduces the data much better than the classical solution W = Wo{M) from Eq. (|32p . The agreement could be 
further improved by taking into account interactions on the interface between the droplet and the fiuid, where the 
fluctuations Arui become non-Gaussian. We will not do this here but instead we observe that subtracting a small 
correction from Wi(M), 

W^2(M) =IFi(M)-2, (51) 

is enough to almost perfectly reproduce the data as shown in Figs. [5] and [71 A physical meaning of this correction 
could be that interactions at the interface droplet-fiuid exert a pressure on the droplet that shifts its boundaries 
towards the centre of mass by one lattice unit on each side of the droplet. 



V. CORRELATED-FLUID PHASE 

In models such as B-in-B [l3| or ZRP [l^ one usually fixes the density p = M /N of particles and takes the 
thermodynamic limit M,N — >■ oo. The condensate emerges in this limit above the critical density pc- The same 
remains true in our model. However, there is another important limit here, namely M^/'^/N = w = const and 
M,N — )■ oo. In this limit, the width W ^ M^^^ of the condensate becomes a finite fraction of the system size N. 
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FIG. 7. Plots of a normalised deviation between theoretical Wth and experimental Wsxp results versus M . Experimental 
results were obtained by MC simulations of the model. Blue symbols show the ratio {Wth — Wexp) /Wexp for classical prediction 
Wth = Wo{M) from Eq. (|32|l . Green and red symbols show the ratio for quantum predictions Wth ~ Wi{M) (calculated 
from Eq. (|50p ) and Wth = W2{M) (Eq. (|51[) ). correspondingly without and with interface corrections. Circles correspond to 
(ci, C2) — (1, 5), squares to (ci, C2) = (2, 5). One can see that the expression W = W2{M) almost ideally reproduces results of 
simulations. 



It turns out that there is a new phase transition as a function of the parameter w: when the width of the condensate 
becomes equal to N, both borders of the condensate merge together. The envelope of the condensate loses its cos'^ 
shape and becomes flat: mean occupation numbers fhi = M/N oc are much larger than 1, and fluctuations which 
are of order -^Z M/N arc not powerful enough to cause {rrii} to drop to mi w 1. Therefore, the condensate no longer 
separates from the background. We shall stress that the existence of this phase is possible only due to periodic 
boundary conditions. If boundary conditions were flxed, i.e., mi = m^v = const, the droplet would not disappear but 
only changed its shape for W > N. 

The correlated-fluid phase is not the same as the weakly-correlated fluid phase below pc- In particular, correlations 
between different m^'s are very strong in this phase. To calculate correlations cov(mj,mfc) = mJrrTk — fhjffik-, let us 
first observe that the partition function ([5|) can be approximated in this phase as 

00 00 
ZcoTT.fiuid{N,M) ^ ... exp 

7711— — 00 TnjY — — CXD 




because the average occupation numbers fhi p and, since we anticipate that yjYai(mi) <ti fhi, we can focus on 
small deviations only. If we now replace the sum by an A^-dimensional integral over mi, . . . , mjv, Eq. ([5^ reduces to 
a Gaussian integral with the constraint on the total number of particles. We can subsequently get rid of the Dirac 
delta by replacing it by 



1 



i5 (x) = lim 

We now define an auxiliary function G{M, N, u) with auxiliary variables u: 



G{M, N, u) 



lim / d^TO 

a-i-O 



1 



27rCT 



exp 



r — — m Am + (b + u) m 
2(7^ 2 



(53) 



(54) 



'^Aij + -^Sij, where Aij denotes a Id discrete Laplacian with periodic boundary 
conditions, and Sa is the Kronecker delta. We have: 



in which bi = Mja^ and A, 
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cov(mj, TTife) = mjmk — mjmk — 
The Gaussian integral in Eq. (j54p can be performed exactly: 



d d 

duj duk 



lnG'(Af, N,u) 



u=0 



G{M, N, u) = lim 



(2^) 



N-l 



(T->o V cr^ det A 



exp 



b^A-'u^ 



-li'^A-h 



(55) 



(56) 
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FIG. 8. Plot of A{k) in the correlated-fluid phase (sohd red hne) calculated from Eq. (|64[) and compared to Monte-Carlo 
simulations for M = 72400, — 80, ci = 1.0, C2 — —0.5 (dashed green line). 



and we obtain that 



var(77ij) — cov{mj ,mj) = lini (A 



cov(mj,mfc) — lini ^) 



' jk 



(57) 
(58) 



The inverse matrix A ^ which appears in these formulas can be calculated using spectral decomposition of the matrix 
A: 



Ajk = ^ Xii>i,jipi,k, 

i 

i 

in which {Aj} and {il^i} are the eigenvalues and the corresponding normaUsed eigenvectors of respectively, 



(59) 
(60) 



N/a^ k = 1 

_ 8ci/p k = 2 

Afc - < 8£i sin2(7r(fc- l)/2Ar) fc = 3,5,7.., 

^sin2(7r(fc-2)/2A^) fc = 4,6,8... 



f 1/ViV 



k = 1 



/JV/2 
sin(-nj(k~2)/N) 

JN/2 



2 

3,5,7, 
4,6,8, 



(61) 



Using the expansion and taking the limit cr ^ wc obtain for large N: 

var(mj) = A//(24ci), 

M ^ 1 (2'Knik-i) 

cov(mj, TOfc) = V — cos - 

4ci7r^ V 

n— 1 ^ 



This means that the correlation function A{k) ~ cov(toi, mfe+i)/var(TOi) behaves as 

All \ 6 1 f 2'!Tnk^ 



n=l 



N 



(62) 
(63) 

(64) 



and does not depend either on ci,C2, or the number of particles. In Fig. [5] we show that A{k) calculated from the 
above equation agrees very well with the result of numerical simulations. 

We shall now discuss the phase transition between the droplet and the correlated fluid phase. For any fixed 
w = M^Z-^/N, there is a critical line in the (ci,C2) phase plane which separates these phases. In the limit of large 
N, M and fixed w, the line can be determined from the condition that W{M) = N (Eq. ([5^ ): 



. ^ 8/3 
C2,crit(Cl) « ITFSTT^ Ci, 



Ar8/3 

where r is the proportionality coefficient r 6.66432 from Eq. 
determined in computer simulations for different M, and rescale ci 



(65) 



This means that if we plot the transition lines 

™ ,2/3 

jV8/3 '^ii ^-ll of them should collapse onto a 



14 



S 2 




0.004 



0.008 0.012 

ciM2/3iV-«/^ 



0.016 



FIG. 9. Rescaled phase line between the droplet phase and the correlated-fluid phase for systems with different number of 
particles M = 4525, 9050, 18100, . . . , 289600, and = 80. Solid black line is our theoretical prediction 
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FIG. 10. Probability of order parameter 7 for droplet-correlated fluid phase transition, A'^ — 80, M — 18100, cl — 0.5, c2 — 1.29. 



single line. We show in Fig. IHlthat such a collapse indeed seems to take place for large system sizes. However, for the 
largest M = 289600 for which we were able to obtain the phase diagram numerically, the data points are still quite 
far from the theoretical line. We believe that this is caused by a very slow convergence towards the asymptotic result 
(j65p due to finite-size corrections which are very strong in the region between the droplet and the fluid. 

The phase transition is of the first order. One can see this by observing that the two phases coexist at the transition 
point with a characteristic binomial structure of the distribution of the order parameter. The double maximum seen 
in Fig. [To] indicates that the system jumps from one phase to another. This is a typical feature of the Ist-order 
transition. 



VI. OTHER PHASES 



We shall now briefly discuss three other phases which appear in our model: localised (condensed) phase, antiferro- 
magnetic fluid, and uncorrelated fluid. For positive ci , C2 , the width W of the droplet decreases with decreasing ci as 
seen from Eq. ([5^ . Finally, at the point of ci = 0, the width formally reaches zero. This means that the condensate 
becomes localised at a single site. For ci ~ 0, the partition function reads 



M 



M 



Z{N,M)^ - E 



mjv = l 



El/3 



M 



(66) 
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and the probability of microstates factorizes over sites, P{{mi}) = /(mi) • • • f{mN), with /(m) = exp(— C2m^/'^). In 
this limit, our model corresponds to the B-in-B/ZRP model with a stretched-exponential weight function /(m) [isj . 
In particular, following [13, Ell, the critical density is given by 

PC = (67) 

with 

oo 

F{z) = z"'exp{-C2m^/^). (68) 

m— 1 

The above series does not admit a closed form, but it can be evaluated numerically for any z and hence the critical 
density (j67p can be computed for any C2 > 0. An important observable in this phase is the distribution of particles p{m) 
- the probability that a randomly chosen node has m particles. This corresponds to the distribution of three- volume 
in CDT. This distribution can be approximated as follows for pc'. 

Pirn) « exp(-C2mi/3)/F(l) + (l/A^)pcond(TO), (69) 

in which the first term corresponds to the critical distribution in the liquid bulk, and Pcond("^) denotes the probability 
of finding m particles in the condensate. We can use the method of Ref. [13, HH to express this probability as follows: 



Pco„d(m) = iV/(m)^ ^JnW) ■ ^^"^ 

Here I{N — 1, m. A/ — m)/Z{N, M) is the probability that the condensate has m or less particles, 



/(7V,m,A/)= Y.--- ^ 



mi — 1 mjv— 1 



/(mi).../(m^). (71) 



Following Ref. |21| . we replace the Delta function by its integral representation, and perform the sum over {mi}. This 
gives 

r°° ds 

I(N -l,m,M -m)^ — exp [TV (ps - ms/A^ + In F(e""))] . (72) 



too 



The integral over ds is dominated by its small-s behaviour. We therefore expand lnF(e '') at s = 0, 

InFie-n . InF(l) - + - ^ + InF(l) - .p. + f! L - P? + ^ V (73) 

^ ' ^ ' F{1) 2 \F{1) F2(l) F(l) J ^ ' 2 y F{1) ^ ' 

and evaluate the resulting Gaussian integral. We obtain that the distribution of mass in the condensate ([70]) is 
approximately Gaussian for m close to N{p — pc): 



Pcondim) cx exp 



2{pc-pi + ^) 



(74) 



This result agrees qualitatively with the simulations, see Fig. [TT] 

When ci < and C2 > 0, the condensate is still localised but the critical density is now zero, i.e., all particles go 
into the condensed phase. This so-called complete (or strong) condensation has its origin in the fact that the radius of 
convergence Zc of Zm{z) from Eq. ^ is becomes zero. This is because Tmn{z) is unbound as either m or n approach 
infinity. The number of particles in the condensate is Af and virtually does not fluctuate. The transition between 
the droplet phase and the localised phase is of second order, because the order parameters are continuous at ci = 0. 

We shall now briefly discuss the antiferromagnetic phase. This phase exists in the region of both coupling constants 
being negative: ci < 0, C2 < 0. The two-point weight f;(m, n) from Eq. ^ has now two positive terms: {rn—n)'^/{m+n) 
which prefers large differences in occupation numbers on neighbouring sites, and m}/^ which prefers large occupations 
but itself does not lead to condensation. In Fig.[T2]we show the correlation function A{k) for this phase. Its oscillatory 
behaviour reflects altered arrangement of occupied/empty sites. Interestingly, the correlation length is quite long. 
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FIG. 11. Probability Pcond(n^) that the locahsed condensate has m particles. Blue symbols: computer simulations for = 
50, M = 2000, ci = 0, C2 = 5. Black line: exact probability distribution obtained from Z{N, M) calculated recursively as in 
Ref. [l^l. Red line: approximate formula (|74|) normalised so that ^^Pcond = 1/A*'. 
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FIG. 12. Gorrelations in the antiferromagnetic phase obtained in MC simulations A'^ = 80, M = 18100, ci 



-0.5, C2 = -0.5. 



which may indicate a possible coupling between two neighbouring occupied sites via a not-completely-cmpty site 
between them. 

Finally, let us consider the uncorrelatcd fluid phase which exists for 01 = 02 = 0. The action 5'[{mi}] equals zero 
and the partition function can be calculated exactly: 



MM/ \ 

m-i—1 mjv — 1 \ ? / 



M - 1 

M -N 



We can now calculate the distribution of particles as follows (cf. Ref. j23|): 

Z{N -1,M -m) _ [N - 1){M - N)\{M - m - l)\ 



p{m) 



1 



Z(7V, M) 



{M -1)\{M -m- N 



exp(-TO/p), 



(75) 



(76) 



where the last formula holds for M N, i.e. for large density p = M/N we typically deal with in this work. The 
distribution of particles (which corresponds to the distribution of three-volume) falls off exponentially with m. 



VII. UNIQUENESS OF g{m,n) 



The choice of the transfer matrix g{m,n) made in Eq. ^ to reproduce the bell-shaped quantum universe is not 
unique. In fact, there is a whole family of functions g{m^n) which lead to the following continuous limit: 



P(mi, . . . , mj^) — > P{m{t)) — exp 



dt ( ci^^^^ +C2m(t)i/3 
m{t) 



and reproduce the shape given by Eq. ([T]). In particular, two other forms of g{m,n), the asymmetric one 

/ \ f {m-nY 1/3 
g[m,n) = exp — ci C2m ' 



(77) 



(78) 
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FIG. 13. Comparison between the average droplet shape rrii (top) and quantum fluctuations y^var(m7) — y m'f — (bottom) 
for asymmetric and symmetric g{m,n), for N = 80, M = 367200, ci = 0.01, C2 — 0.59. Black curves: MC simulations of the 
full CDT model (courtesy of A. Gorlich) for T — 80 time slices and the total volume 14 — 367200 equivalent to the number of 
sites and particles in our simulations. Left: symmetric g{m,n) from Eq. (|79|l . The same result is obtained for the symmetric 
Eq. ((4| studied in previous sections. Small asymmetry in ^/^mr(Jni) is caused by statistical fluctuations. Right: asymmetric 
g{m,n) from Eq. dtSJ. 



and the symmetric one with the geometric mean ^Jmn rather than the arithmetic mean (m^Ti)l^ in the denominator, 

5(m, n)=exp -ci ^= C2 , (79) 

\ \Jran I ) 

have the same asymptotic behaviour as Eq. Our simulations show (see Fig. [T5)) that the shape of the droplet is 
reproduced well by all three forms of f/(m, n) in the large- A^, limit. However, the shape is slightly asymmetric in 
the case of Eq. (j78p . whereas it is perfectly symmetric for symmetric forms of g(m^n) as those given in Eqs. (U) or 
([7^ . However, the data from the full CDT model are perfectly symmetric (excluding small statistical fluctuations). 
We thus conclude that the asymmetry is of finite-size origin and that the effective transfer matrix in CDT has to be 
symmetric as in Eqs. (|4|) or ((79|) . Interestingly, although Eq. ((79)) leads to exactly the same envelope (HJ in the droplet 
phase as Eq. it does not permit the existence of the antiferromagnetic phase. Indeed, the corresponding action 
in the antiferromagnetic phase, 

5a„tiferr. ~ Ic^R-^ ' ' ^ + C-.K'' I ^ ' , (80) 

is bigger than the corresponding action in the localised phase, 

^locali.cd « 2ciAf3/2 + C2MI/3, (81) 

for ci < and for any C2, and therefore antiferromagnetic states are disfavoured in this case. In other words, the 
localised phase extends to all C2 (positive and negative) in the phase plane (compare with Fig. 2]) for the model with 
the transfer matrix given by Eq. (j79|) . We see that the existence of the antiferromagnetic phase depends on the 
behaviour of the kernel g{m, n) for small values of the arguments. 



VIII. CONCLUSIONS 



In this paper we have analysed a simple model of particles residing on sites of a Id lattice, in which the probability of 
microstate ^ equals , where 5* corresponds to the effective action ^ of the CDT model. We have shown that our 
model reproduces not only the average shape of the droplet - the macroscopic universe of CDT ~ but also quantum 
fluctuations around it. We have calculated the extension of this droplet and shown that the quantum universe is 
bigger than classical de-Sitter solution. 
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The droplet phase is one of five different phases whieh exists in our model. Two of these phases, localised condensate 
and uncorrelated fluid, can be identified as phases "B" and "A" of CDT. In each of these phases, we have calculated 
the distribution of particles p{m) which corresponds to the distribution of three-volume in CDT. By measuring this 
distribution in the original CDT model and comparing it to our predictions one could validate our hypothesis that 
all phases can be described by the same effective action. 

Furthermore, we have predicted the existence of at least one more phase - the correlated fluid phase. This phase, 
although yet unobserved, must surely exists in CDT as a simple consequence of periodic boundary conditions ensured 
by the global topology of CDT. We have calculated two observables: p{m) and the correlation function A{k), which 
can be easily measured in CDT. The agreement with our predictions would provide further evidence for the effective 
action ([2])- 

Lastly, we have suggested that, depending on the behaviour of the action for small three- volumes, the fifth, anti- 
ferromagnetic phase can exist. 

Our predictions can be tested in the CDT model, even without the knowledge of the mapping between the effective 
coupling constants ci, C2 and the parameters in the Einstein- Hilbert action of CDT. In particular, the values of ci, C2 
can be determined by fitting Eq. pop to the data from computer simulations in the macroscopic-universe phase, 
calculating Ci/c2 from W, and resolving for ci,C2 using the equation for the background density pc- Then, the 
correlated-fluid phase can be reached by increasing M. In other phases, equations derived in this paper for some 
quantities can be used to determine ci,C2. These values in turn can be applied to calculate other quantities and 
compare them to those estimated in full CDT simulations. 
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APPENDIX - NUMERICAL SIMULATIONS 



Our model can be simulated using standard Monte Carlo techniques. We start each simulation from some initial, 
random configuration of particles and construct a Markov chain in the space of configurations by moving particles 
between sites with probability depending on the current configuration. More specifically, we construct a new config- 
uration B = {nil, . . . , — 1, . . . , nij -|- 1, . . . , m^} from the old one A = {mi, . . . , m^, . . . , nij, . . . , mjv} by picking 
two random sites i and j with > 1, and moving one particle from site i to site j with probability given by the 
Metropolis formula 



P{A ^B)= min (l, ^\ = min (l, si^^-^^^'; ' ; m.^i).9(m,_i, + l)g(m, + 1, m.+Q | ^g^) 



if i, j are not nearest neighbours, and with probability 

uf A . u\ • /i "T-i - l)g('-»» - 1, -I- l)g(m,+i + 1, m,+2) ) 

P{A -> B) mm <^ 1, — r }■ forj = i + 1, (83) 

1^ g(^ni^l,mi)g{mi,rni+l)g{■mi+l,m^+2) J 

uf A . u\ • /i .g(™^~2, + l)g{mi-i + l,mi- l)g(m, - 1, m,+i) ) . . 

P(A B) = mm <^ 1, ; — — r > forj = i - 1, (84) 

|_ g{mi-2,mi-l)g{rni--^,rni)g[m^,m^+l) J 



if they arc neighbours, i.e., if \i — j\ = 1. Such form of the acceptance probability guarantees that the probability of 
microstate P{{mi}) will be given by Eq. ([3]). It is convenient to introduce the following notation: 

g(m,n- 1) g(m- l,n) g{m-l,n+l) g{m + l,n-l) 

a(m,n) = — , fi(m,n) ^ — , 7(m,n) = , d(m,n) = r . (85) 

g{m,n) g[m,n) g(m,n) g[m,n) 



19 



Then, the acceptance probabilities can be rewritten as: 

PM i?) = mm <^ 1, — — — ; ->, fork-j >1, (86) 

^ ^ i a(mj_i,mj + l)/3(TOj + l,mj+i)J i i ' \ ' 

P{A B) = mm <^ 1, — r \ , for = i + 1, (87 

[ /3(mj+i + l,mi+2) J 

P(A^P)=mmn, — }, forj=i-l. (88) 

(89) 

In our simulations, we calculate and store the values of a{m, n), l3(rn, n), "f(m, n), S{m, n) for m, n = 1, TOmax, with 
some m„iax < M. This allows us to use Eqs. and to avoid time-consuming computations of the ratios of 

g{m,n) in Eqs. if only the number of particles at sites i,j does not exceed TOmax- Otherwise, we calculate 

the acceptance probability directly from Eqs. (|82|) - (|84| . The value of mmax - typically a few thousands - is chosen as 
big as possible given available computer memory. To reduce the autocorrelation time, measurements are made every 
M moves. 

All measurements of the average shape of the condensate and fluctuations around it are performed by shifting the 
condensate for each sample to a common centre of mass at site i = N/2. In order to account for periodic boundary 
conditions, the centre of mass is found in a 2d plane, assuming that the sites reside on a circle in this plane, and then 
the coordinates {x, y) of that point are mapped to the index i of a site closest to the centre of mass. We have checked 
that other procedures of finding the centre lead to very similar results. 
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